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ABSTRACT 

The standard equilibrium for radiation-dominated accretion disks has long been known to be viscously, 
thermally, and convectively unstable, but the nonlinear development of these instabilities — hence 
the actual state of such disks — has not yet been identified. By performing local two-dimensional 
hydrodynamic simulations of disks, we demonstrate that convective motions can release heat sufficiently 
rapidly as to substantially alter the vertical structure of the disk. If the dissipation rate within a vertical 
column is proportional to its mass, the disk settles into a new configuration thinner by a factor of 
two than the standard radiation-supported equilibrium. If, on the other hand, the vertically-integrated 
dissipation rate is proportional to the vertically-integrated total pressure, the disk is subject to the 
well-known thermal instability. Convection, however, biases the development of this instability toward 
collapse. The end result of such a collapse is a gas pressure-dominated equilibrium at the original column 
density. 

Subject headings: accretion, accretion disks, convection 



1. INTRODUCTION 

More than twenty-five years ago, Shakura & Sunyaev 
(1973, hereafter SS) showed that the inner portions of ac- 
cretion disks with luminosities near Eddington are likely to 
be dominated by radiation pressure. Because the vertical 
component of gravity increases oc z near the disk mid- 
plane, true vertical hydrodynamic equilibrium could only 
be achieved if the heating rate is constant with height. If 
the physical dissipation rate scales locally with the density, 
then hydrostatic and radiative equilibrium require the gas 
density to be constant from the disk midplane to the disk 
surface. 

Unfortunately, this equilibrium is subject to numerous 
instabilities. Lightman & Eardley (1974) pointed out that 
if the viscous stress is proportional to the radiation pres- 
sure, perturbations to the surface density grow on the 
(comparatively long) viscous inflow timescale. Shakura 
& Sunyaev (1976) then observed that in these same condi- 
tions the thermal content of the disk is likewise unstable, 
with a growth time comparable to the (shorter) thermal 
timescale. Bisnovatyi-Kogan & Blinnikov (1977, hereafter 
BKB) noticed that if the radiation is locked to the gas 
even on short length scales (i.e., if, for the purpose of dy- 
namics, the optical depth is treated as effectively infinite), 
such disks should be convectively unstable, for the spe- 
cific entropy decreases upward (the radiation pressure de- 
creases upward while the density is constant); the linear 
growth rate for convective "bubbles" was worked out by 
Lominadze & Chagelishvili (1984). This work has been re- 
cently extended by Pietrini & Krolik (2000), who derived 
a hydrodynamic WKB dispersion relation in the presence 
of rotation and including the effects of finite optical depth. 
In other recent work, Gammie (1998) demonstrated that 



a magnetic field in radiation-supported disks can catalyze 
a short- wavelength (kH ^> 1, where H is the disk scale 
height) overstable wave mode. Blaes & Socrates (2000) 
found the dispersion relation for radiation MHD modes 
within radiation-dominated disks. 

With so many instabilities potentially operating, one 
must wonder what is the real state of these disks. Numer- 
ous suggestions exist in the literature. BKB argued that 
the disk structure could be determined by imposing the 
conditions of constant specific entropy (a state achieved as 
a result of convection) and hydrostatic equilibrium. Liang 
(1977) defined a new equilibrium by parameterizing the 
efficiency of convective heat transport. Shakura, Sunyaev, 
& Zilitinkevich (1978) found an analytic solution for the 
vertical structure based on these same assumptions, but in- 
cluding support from turbulent motions in the hydrostatic 
equilibrium. They were able to find an analytic solution to 
the equations granted the assumption that the rms Mach 
number of the turbulence took a special value. Robertson 
and Tayler (1981) argued that convection does not quench 
the thermal instability, making rough estimates for the ef- 
fect of convection on transport. Others (Cannizzo 1992; 
Milsom, Chen & Taam 1994; Rozanska et al. 1999) have 
proposed models based on the "mixing-length" prescrip- 
tion. Eggum et al. (1987), Milsom & Taam (1997), and 
Fujita & Okuda (1998) used 2-D radiation hydrodynamic 
simulations to study various aspects of these disks. 

Unfortunately, each of these previous attempts to un- 
derstand the structure of radiation-dominated disks has 
been lacking in one or more crucial respects. None of the 
analytic efforts actually solved the force equation with- 
out making some assumption about the character of the 
answer, while none of the numerical simulations had the 
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resolution to actually determine the internal vertical struc- 
ture (these simulations were primarily aimed at exploring 
issues involving the global behavior of disks rather than 
the particulars of vertical structure). 

The object of the work presented here is to use radi- 
ation hydrodynamic simulations to focus on the vertical 
equilibrium of radiation-dominated disks so that one may 
actually determine the structure of bright accretion disks 
in the range of radii where they release most of their en- 
ergy, their innermost rings. In addition to the intrinsic 
interest of this effort for deepening our understanding of 
accretion dynamics, solving this question is a prerequisite 
for any effort to predict the spectrum of radiation emerg- 
ing from these disks — the nature of the disk atmosphere, 
and therefore the character of any features imprinted on 
the spectrum in the disk photosphere, depends crucially 
on the vertical distribution of gas density and heat depo- 
sition. Though the viscous mechanism operating in accre- 
tion flows is likely due to magnetic stresses, we parame- 
terize the effects of viscosity so as to use a simpler hydro- 
dynamic code. 

2. SIMULATION METHODS 

The tool we employ is the Zeus simulation code (Stone 
& Norman 1992) with its flux-limited radiation diffusion 
module (Turner & Stone 2001). This code solves five cou- 
pled partial differential equations on a fixed Eulerian grid: 
the mass-continuity equation, the Navier-Stokes equation, 
the energy conservation equations for both the gas and 
the radiation, and the radiation momentum conservation 
equation. Defining the radiation quantities in a cylindrical 
coordinate frame comoving with the fluid, and retaining 
only those terms first order in v/c, these equations are: 
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Here the convective derivative D/Dt = d/dt + v • V. The 
unit vectors in the radial, azimuthal, and vertical direc- 
tions are r, <fi, and z, respectively. The quantities p, e, 
v, and p are the gas mass density, energy density, veloc- 
ity, and pressure respectively, while E, F, and P are the 
radiation energy density, momentum density or flux, and 
pressure tensor, respectively. Only absorptive opacity is 
included in the Planck mean opacity Kp and the energy 
mean opacity ke, but scattering is added to absorption in 
the flux-mean total opacity k. Energy injection is permit- 
ted via the function rj (see later discussion). The orbital 
frequency, fl — yjGM/r^, is evaluated at the central ra- 
dius r Q . The disk is assumed to be azimuthally symmetric, 



and shearing-box coordinates (the Hill potential) arc used 
to describe the gravity (the ^-dependent terms in equa- 
tion 2 all arise from this approximation to the potential in 
rotating coordinates). 

The mass-continuity and Navier-Stokes equations are 
advanced in time using operator-splitting. The radiation 
transport problem is solved implicitly using the approx- 
imation of flux- limited diffusion (Turner & Stone 2001). 
The code assumes a gas equation of state with e, the gas 
energy density, evolved adiabatically, but including ab- 
sorption and emission of radiation. The gas pressure is 
defined as p gas = (7 — l)e, with 7 = 5/3. 

In all our simulations, the problem area was a radial 
segment of a geometrically-thin, optically-thick, radiation- 
dominated accretion disk. In terms of the disk height, 
H = F kf/{cQ 2 ) as predicted by the SS equilibrium, 
where F Q = 3Q 2 M/(8n) is the radiation flux at the top 
of the disk, c is the speed of light, we simulated the dy- 
namics in a region stretching in the vertical direction from 
z = —2H to z = 2H, and in the (cylindrical) radial di- 
rection from r — r Q — 2H to r = r Q + 2H. Two values 
of central radius r a were chosen, 100r 9 and 200r g , where 
r g = GM/c 2 . The black hole mass was 10 8 M Q , and we set 
the accretion rate at M = (3 — 10) x 10 25 gm s _1 , appro- 
priate for bright active galactic nuclei shining at (20-60)% 
of the Eddington limit. 

The initial condition in every case was a slightly- 
modified SS equilibrium. Because that equilibrium be- 
comes ill-defined at the disk surface (where it predicts 
that the gas density falls discontinuously to zero) , we con- 
structed a more complete version of the same equilibrium 
allowing for the small amount of gas pressure support that 
exists at the top. In this initial equilibrium (but not in our 
simulations), we assume that T gas — T rac i locally through- 
out the disk. Fixing the flux, F Q , column density, E, and 
orbital frequency, O, we used a shooting method to solve 
simultaneously the hydrostatic equilibrium equation, 
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the radiative equilibrium equation, 
dF _ 2F oP 

~dz ~ E ' 
and the radiation diffusion equation, 

p _ C dp rad 

[up + ot.o) dz 

We introduced an scattering opacity floor, a Q , to prevent 
the opacity from ever becoming small enough to under- 
cut the validity of the diffusion approximation (see further 
discussion below) . The gas density in this solution is very 
nearly constant from z — almost to z = H; starting 
from just below z = H, it falls steeply, but not discon- 
tinuously, as z increases. The radial and vertical velocity 
components were set to zero initially, while the azimuthal 
velocity components were set to the shearing sheet value, 
-1.5fi(r-r„). 

The column density of the accretion disk at the start 
of the simulation is chosen to be the value in the the SS 
equilibrium. Ignoring relativistic correction factors and 
the correction factor accounting for the outward angular 
momentum flux, it is 

E - J-rh^x 3 / 2 , (9) 
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where m is the accretion rate in Eddington units (for unit 
efficiency), a is the ratio of the viscous stress to the pres- 
sure, £ has units gm cm~ 2 , and x — rc 2 /GM. In all 
but one special case, we fixed a at 0.01 for computing the 
surface mass density at the start of the simulation. 

At the top and bottom edges [z = ±2iT), the boundary 
condition was chosen to be outflow; that is, fluid quanti- 
ties in the cells adjacent to the problem area were set equal 
to the boundary values. In addition, anywhere we set an 
outflow boundary condition we required the velocity to be 
outwards; if not, it was set to zero. In the radial direction, 
periodic boundary conditions were assumed in all quanti- 
ties except v,/), for which periodicity was enforced for the 

quantity Ity - V K eplerian- 

The scattering opacity was taken to be the electron scat- 
tering value K es = 0.4cm 2 gm -1 , with certain exceptions. 
When regions of the simulation zone were optically thin, 
we found that the diffusion routine took many steps to 
converge, greatly slowing progress of the code. In ad- 
dition, the outflow boundary condition for optically-thin 
flux-limited diffusion was impossible to implement due to 
a steep density gradient and uniform radiation field near 
the boundary where the disk has very low density but is 
nonetheless gas pressure supported. To avoid these dif- 
ficulties associated with optically-thin radiation transfer, 
an opacity floor was chosen such that the scattering op- 
tical depth across each cell was at least unity (a was 
2 x 10 _12 cm _1 for simulation 1 and 2 x 10~ 13 cm~' for sim- 
ulation 1). This adjustment of the opacity affected about 
half of the simulation region, but a much smaller fraction 
of the total mass. It caused the radiative flux through the 
upper and lower boundaries to be carried advectively with 
the fluid rather than by diffusion relative to the fluid; this 
forced outflow through the outer boundaries, maintaining 
consistency with our outflow boundary condition. Over a 
characteristic thermal time for the Shakura-Sunyaev equi- 
librium 

tthermai = F' 1 / (E + e)dz = HkV/c ~ 4/(fia), (10) 
Jo 

only a small fraction (<; 1%) of the mass within the 
simulation region was lost due to outflow. The absorp- 
tion opacity was assumed to be purely bremsstrahlung: 
KabsP = 10 52 p n / 2 e~ 7 / 2 cm -1 , where p is the mass density 
and e is the gas energy density. 

Without creation of new photons, flux lost out of the top 
of the disk would ultimately deplete the disk of radiation. 
In a real disk, new photons are created in a way that de- 
pends on the local gas density, temperature, and magnetic 
field. Here we wish merely to model the radiation process 
in a phenomenological fashion. To this end, we simply call 
the local radiation rate rj (units of erg cm. s^ 1 ) and de- 
fine it by either of two prescriptions: in one case, ij oc p, 
while in the other r/ is proportional to the total pressure 
averaged over the simulation region, but locally propor- 
tional to the gas density within any given cell. The for- 
mer choice was meant to mimic the assumption underlying 
the SS vertical equilibrium. The latter choice is based, of 
course, on the thought that the stress is proportional to 
the total pressure. 

Ultimately, the energy for these new photons comes from 
dissipation of gas motions (and, in a real disk, resistive 
dissipation of magnetic field energy). However, there is 



no clear-cut way to connect photon creation to local vis- 
cosity. For example, if, as seems likely, most of the angu- 
lar momentum transport in disks is due not to something 
that behaves like viscosity, but rather to MHD turbulence 
(Balbus & Hawley 1998), there is no simple, local connec- 
tion between stress and heating, much less between stress 
and radiation. Consequently, we do not place any viscous 
counterpart to the radiation creation term rj in the gas 
momentum and energy equations. 

The final item to be noted in specifying the simulations 
is the resolution. The results we present were done with a 
128 x 128 grid, but we reran part of simulation 1 with reso- 
lutions of 64 x 64 and 256 x 256 to check for convergence. In 
those cases in which the disk collapsed as a result of ther- 
mal instability, we stopped the simulation when nearly all 
the disk mass was contained within then central 1/4 of 
the vertical coordinate. At that point, the central half (in 
both radial and vertical directions) was rebinncd to twice 
the resolution, and the simulation was then restarted. 

The simulations and their characteristics are listed in 
table ^. The three parameters varied were the radius, 
M, and resolution; the radius and M change the optical 
depth as described by equation [| and also change the ratio 
of Prad/Pgas- Each simulation consisted of two successive 
phases: 

A) The disk was first allowed to come to a statistically 
steady equilibrium modified by convection with the inte- 
grated heating rate fixed globally, but scaled locally with 
the density, r\ = F p/H. This readjustment took several 
thermal times. 

B) The heating rate was then allowed to vary in pro- 
portion to the average total pressure ( "a" prescription) to 
search for possible thermal instability. In this case 
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where (p ) is the volume averaged pressure at the end of 
Phase A. 

Table 1: Simulation Parameters 



N 


M(g/s) 


r/r g 


Res. 


4H (cm) 


4>rad 


<t>gas 


a<, 


1 


3.E25 


100 


128^ 


2.4E13 


95 


0.50 


2.E-12 


2 


1.E26 


200 


128 2 


1.6E14 


160 


0.95 


2.E-13 



3. READJUSTMENT BY CONVECTION 

In this section we discuss the results from Phase A of the 
simulations. 



3.1. Approximate Analytic Description of the 
Time-Steady State 

To provide a context for the simulation results, we begin by 
discussing how one might find approximate analytic solutions 
for the vertical structure of radiation-dominated disks. Because 
the inflow time is very long compared to the dynamical time, 
tdyn = 2ttQ~ 1 , one might expect the disks to be in hydrostatic 
equilibrium. In the absence of convection, one would also ex- 
pect the heat flux to be carried by photon diffusion. 

On the other hand, when convection is active, one might 
estimate the heat flux in terms of mixing- length theory (e.g., 
as described in Clayton 1968). If one could define an effective 
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mixing-length the fraction of the heat carried by convection 
is 

F conv /F a = 13(z//i)" 1/2 (///i) 2 (a/0.01)(/iAVlnT) 3/2 , (12) 

where A V In T is the difference between the true logarithmic 
temperature gradient and the adiabatic temperature gradient. 
To the degree that convection is efficient, entropy gradients are 
erased, so that the disk becomes nearly isentropic. 

Thus, one way to approximately determine the structure of 
these disks is to require all three of these conditions: hydro- 
static balance, isentropy, and photon diffusive equilibrium (or 
at least that photon diffusion carry a specified fraction of the 
heat). However, this is impossible because any two of these 
three conditions suffice to determine the structure of the disk. 
For example, SS applied the two conditions of hydrostatic equi- 
librium and photon diffusive equilibrium, without reference to 
the condition of isentropy (in fact, their equilibrium, as we have 
already remarked, has a strongly unstable entropy gradient). 
In this solution, the density is constant with altitude up to H, 
and drops sharply to zero for z > H. Thus, the mass- weighted 
average height of the disk is (z) = J pzdzj J Q pdz ~ H/2. 

Similarly, BKB assumed hydrostatic balance and isentropy, 
but made no assumption with regard to radiative balance. If 
the turbulence is subsonic, then this approximation is appro- 
priate since the turbulent pressure will be much smaller than 
the radiation pressure. Then, the diffusive flux is simply 



F dif = F z/H. 
The density distribution they found is: 
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where p c is the central density, given by 
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The mass-weighted average height of this disk is (2) = 
35///128, about half that of the SS solution. The specific en- 
tropy, s = |a 1/ ^ 4 _E 3/ ' 4 /p, is found to be 
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In this case, if the dissipation is proportional to density, then 
the convective flux is 
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where £ = z/H. Thus, for small £, the ratio of the convective 
flux to the diffusive flux is simply 19/16; i.e. the convective flux 
is always comparable to the diffusive flux. At £ = 1, F co „ = 
so all the radiation escapes diffusively. The volume-integrated 
radiation energy in the disk is reduced by a factor of 3 with re- 
spect to the SS solution, so the thermal timescale is decreased 
by the same factor. 

Alternatively, one could also insist that the specific entropy 
is constant and that the disk be in radiative balance but pay 
no attention to hydrostatic equilibrium (this assumes that the 
convective flux is small compared to the diffusive flux). In that 
case, if the dissipation is proportional to density, then the den- 
sity profile would be 
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where cn is a Jacobi elliptic function and p c is the central den 
sity, given by 
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The entropy turns out to be identical to equation 
[lil except for a different numerical factor in front: 
4(4/3) 1/4 r(3/4)/r(l/4)/0F, which is 1.3 times greater than 
in the BKB solution. The average radiation energy density 
in this solution is 4 times smaller than in the SS solution. 
A disk obeying these assumptions would be out of hydro- 
static balance by p c S7 2 i/(v / 2dn(ujl/2)sn(u|l/2) — z/H), where 
u = \/2pc^/E, and dn, sn are Jacobi elliptic functions. For 
hydrostatic equilibrium to exist, turbulence must provide this 
deficit. The quantity in parentheses can be approximated by 
~ 0.4shi7r(z/i/) 4 / 5 , so the turbulent velocity must be of order 
the sound speed, ~ HQ. 

Thus, at most two of these three plausible conditions can 
be satisfied exactly. On the other hand, it is possible to ap- 
proximately satisfy all three if small departures are allowed for 
each one. Shakura et al. (1978) tried to modify both the hy- 
drostatic balance and energy transport equations to allow for 
small departures (specifically, the kinetic energy contribution 
in the hydrostatic equilibrium equation and convective heat 
transport in the energy equation). However, they were able 
to find an analytic solution only for the special case in which 
the rms Mach number A4 = 2/3. When that is the case, the 
density distribution is 



p = pc 



(20) 



where p c = 3E/(4_H'). Unfortunately, there is no particular rea- 
son to expect that M takes this value (in fact, the simulations 
we performed indicated that it is smaller by an order of mag- 
nitude). Consequently, it is necessary to solve the equations 
of motion numerically in order to find the true compromise 
between these three conditions that is found in real disks. 

3.2. Dynamical readjustment due to convection 

We are now ready to report the results of our simulations. 

Within the body of the disk, the timescale for thermal equi- 
libration of the radiation and the gas is the shortest timescale; 
it is a fraction < 1/20 of a dynamical timescale. Consequently, 
deep inside the disk the gas and radiation temperatures quickly 
become equal. Near the surface, the timescale for radiation to 
diffuse outwards is comparable to or shorter than the thermal 
equilibration timescale, so in those regions the gas temperature 
is anywhere from 10 to 10 3 times greater than the radiation 
temperature. 

Though our initial condition satisfies hydrostatic equilib- 
rium, it is convectively unstable. Pietrini & Krolik (2000) 
showed that in the WKB approximation (i.e., kH 3> 1) the 
most rapidly growing modes in the linear regime are those 
whose wave-vectors are nearly horizontal, and that the growth 
rate for these modes is 
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almost independent of wavelength when kE 3> 1 (note the ty- 
pographical error of a factor of 2 in equation 20 of Pietrini 
& Krolik). That is, they predicted that convective motions 
in this equilibrium will grow most rapidly near the surface of 
the disk, and that the basic growth timescale is the dynamical 
time. Our simulations confirmed this prediction quantitatively 
for z/H < 1/2; for example, at 0.4H the analytic growth time 
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is 9xl0 5 s, while the measured growth time in simulation 1 is 
10 6 s. However, near the disk surface, the WKB approxima- 
tion breaks down because there are sharp gradients in density 
and pressure. As a result, the numerical growth rate we find is 
slower than the analytic value by a factor of a few near z ~ H. 

In the non-linear regime, growth is most rapid for longer 
wavelengths, much as in the case of the Rayleigh- Taylor insta- 
bility (Garabedian 1957). Horizontal density modulations near 
the top of the simulated disks with wavelengths of order H /10 
grow within a few dynamical times of the start of the simulation 
(Figure Overdense blobs fall, become Kelvin-Helmholtz un- 
stable, and form an inverse-mushroom shape. The disk rapidly 
becomes turbulent, mixing high and low entropy regions, and 
then smooths out a bit, while steady convective motions con- 
tinue. 




Fig. 1. — The density distribution in simulation 1 (at 256 2 res- 
olution) at times t/t dyn = 0, 1, 2, and 3. The grey scale is linear 
from p = (black) to 6 X 10 — 9 gm cm -3 (white). The horizontal 
axis is radial distance from the central radius of the simulation, the 
vertical axis is altitude from the disk midplane; both are in cm. 

After ~ 90 dynamical times, the disk achieves an approxi- 
mate statistical steady-state. Convection continues to occur, 
creating small-amplitude fluctuations, but average properties 
remain very nearly constant. The outgoing flux at the top of 
the box fluctuates about the mean heating rate, evidence that 
the simulation has reached a steady state. Several stable con- 
vective cells are formed, in which low entropy sinking regions 
of size ~ H are surrounded by upwelling plumes with twice the 
entropy (Figure [| to be precise, at any given height within the 
body of the disk, the specific entropy in the plumes is roughly 
twice the lowest specific entropy elsewhere in the disk at that 
height). The plumes are not completely resolved since at their 
narrowest in the midplane, they are only spanned by 2 pix- 
els. If the heating rate is proportional to the local density, 
mean conditions in the disk remain constant for the duration 
of Phase A of the simulation, 11 initial thermal times, or 180 
dynamical times. In contrast to the "square-wave" density pro- 



file of the initial equilibrium, the mean density profile exhibits 
a smooth fall from the midplane out to the top. The mass- 
weighted average height of the disk changes by a factor of ~ 2 
as convection carries more flux outwards, reducing the radia- 
tion pressure, causing the disk to collapse. The reduced height 
further reduces the thermal time, causing more flux to escape, 
and further collapse until radiation pressure supports the disk 
again. At the end of the readjustment process, the volume- 
integrated radiation energy is reduced relative to what it was 
in the initial equilibrium by a factor of 4.4 and 4.6 in simula- 
tions 1 and 2, respectively. This factor is slightly larger than 
the factor of 3 predicted by the analytic BKB scaling, but is 
close to the factor of 4 predicted by equation [l8|. We reran sim- 
ulation 1 at lower (64 2 ) and higher resolution (256 2 ) to check 
for convergence. We found that the radiation energy density 
in the 64 2 simulation disagreed by 20% at the end of phase A, 
demonstrating that at this resolution the simulation was not 
converged. At higher resolution, 256 2 , we were only able to 
run the simulation for about a thermal time since the run time 
scales as N 4 , but we found that the difference between the to- 
tal radiation energy density of the 128 2 and 256 2 simulations 
had a standard deviation of only 0.6%, indicating that the 128 2 
simulation was indeed converged. We utilize this higher reso- 
lution run when studying the initial convective collapse (e.g. 
Figure [j]). 
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Fig. 2. — The specific radiation entropy distribution in simulation 
1 (resolution 128 ) in units of cm 2 s — 2 K — 1 at t = 97t dyn . Axes 
are as in Fig. 1; regions above the maximum of the color scale are 
white. 

The mean time-steady density profile of simulation 1, Phase 
A, is shown in Figure H, where it is contrasted with the var- 
ious approximate analytic solutions discussed in the previous 
subsection. This distribution lies somewhere between the BKB 
and equation IIm solutions. These analytic solutions break down 
near the surface where, in actuality, gas pressure support be- 
comes important and LTE no longer holds. We do not find any 
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dependence of the inner disk structure on the opacity floor (we 
increased it by a factor of 2 and found the disk structure was 
identical). The Shakura et al. (1978) solution, with column 
density equal to the simulation result, deviates significantly 
from the simulation mean. 
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ployed in Zeus. In simulation 1, Phase A, dissipation due to 
artificial viscosity comprises 9% of the total heating rate within 
the box. In the immediate vicinity of a shock, the dissipation 
can be significant, but averaged over larger volumes and times, 
it does not qualitatively alter the disk's thermal properties. 
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Fig. 3. — Density distribution at t = (flat-topped solid line) 
and t = 90t t {y rl in simulation 1 (bell-shaped solid line) compared to 
the prediction of BKB (dotted line), equation 18 (dashed line), and 
Shakura et al. (1978) (dashed-dot line). 

3.3. Heat Transport 

The most dramatic feature in Figure ^| is the sharp contrast 
in specific entropy between rising and sinking regions. Nonethe- 
less, as one might expect from the nonlinear development of an 
instability that feeds on vertical entropy gradients, convection 
does an excellent job of smoothing out the radially- averaged 
specific entropy profile (Figure ^|). Within a few dynamical 
times, the radially- averaged specific entropy in the disk be- 
comes very nearly constant for z < 0.8H; this (on-average) 
isentropic region contains 98% of the disk mass. Convection 
is so efficient that this is accomplished with relatively small 
amplitude motions: the rms Mach number in simulation 1 
is A4 ~ 0.1, vindicating the quasi-hydrostatic approximation 
made in the BKB solution. 

Although the radially-averaged density distribution in the 
disk is close to that predicted by the two isentropic models (see 
Fig. ^), the actual value of the mean specific entropy in the disk 
is about 25% smaller than the analytic prediction (equation 
|l6| ). This departure may be due to the fact that the entropy 
is set by the requirement of radiative diffusion matching the 
total flux at the top of the disk. In the simulations, the density 
gradient near the disk surface is somewhat different than that 
predicted by the analytic models because turbulent support 
becomes comparable to radiation support (A4 ~ 1). As a re- 
sult, the entropy necessary to maintain the flux can be slightly 
different. 

This strong turbulence near the disk surface finds its ulti- 
mate origin in the reduction in optical depth across convective 
cells as the density falls sharply near the disk surface. When 
photons can diffuse across a cell in an eddy turnover time, the 
radiation does not effectively couple to the fluid, and the sound 
speed relevant to fluctuations on those lengthscales drops to the 
(much lower) gas sound speed. Because photon diffusion pre- 
vents the radiation from being compressed along with the gas, 
the gas's compressibility increases dramatically. That permits 
large local density enhancements, which lead in turn to rapid 
cell-sinking and high speed motions. 

Compressive motions call into play the artificial viscosity em- 



FlG. 4. — Horizontally-averaged entropy as a function of height at 
t=0 (solid line), and at t = 95t dyn (dashed line) from simulation 1. 

The mean diffusive flux, Fdif = {cd z E /(3k)) and the con- 
vective flux, F con = (v s {E + e) + f* dz(Vv : P + pV ■ v)>, 
averaged radially and in time, are shown in [B|. The convective 
flux is comparable to the diffusive flux near the midplane of 
the disk, as predicted by the BKB solution. It is, of course, the 
additional heat loss due to convection that reduces the thermal 
time so much. In line with intuitive expectation, the diffusive 
flux rises steadily with height within the disk until it carries 
almost all the heat at the disk surface. The integrated heating 
rate and the total flux are nearly equal at the top of the box, 
indicating that the disk is in quasi-equilibrium. 

With these observations in hand, we can evaluate the pos- 
sible relevance of estimates made through mixing-length the- 
ory. On the one hand, the fact that the mean specific entropy 
is very nearly constant as a function of altitude is consistent 
with the usual expectation that convection very nearly erases 
any difference between the actual temperature gradient and the 
adiabatic temperature gradient. Moreover, our finding that the 
convective and diffusive heat fluxes are similar deep inside the 
disk could in principle allow us to turn around equation |l^ 
and evaluate an effective mixing length. However, there are 
both technical and conceptual problems preventing this. The 
technical problem is that the limited numerical accuracy with 
which we can measure the mean temperature gradient will lead 
to a highly uncertain estimate of the difference between it and 
the very similar adiabatic gradient. The conceptual problem 
is that the large instantaneous contrasts in specific entropy at 
fixed altitude demonstrate that the mixing length picture has 
questionable validity here because the cell sizes are comparable 
to the entire thickness of the disk. 

Above the surface of the disk, the opacity floor that we im- 
pose creates a numerical artifact in the character of heat trans- 
port. With the true opacity, that region would be optically 
thin and the heat would be carried almost exclusively by free- 
streaming radiation; the opacity floor causes the heat to be 
carried advectively. We believe that this has little consequence 
for dynamics in the disk body because the region in which the 
opacity floor is implemented contains only ~ 1% of the disk 
mass, and because a simulation with an opacity floor twice as 
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large gave similar results. 
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Fig. 5. — Horizontally-averaged diffusive energy flux (dotted line) 
compared to the convective flux (solid line), total flux (dashed line), 
and integrated heating rate (dot-dash line) from t/t dyn = 90 to 
95 in simulation 1. The dominance of convective heat flux for 
\z\ > 4 X 10 13 cm is an artifact of our opacity floor; in real disks the 
flux at high altitude is free-streaming. 

4. EVOLUTION ON THE THERMAL TIMESCALE 

4.1. Thermal Instability 

One of the important questions motivating this study is 
whether the thermal instability predicted on the basis of the q- 
model actually occurs in radiation-dominated disks. We began 
this part of our study by verifying that our simulation code 
reproduced this instability in 1-D when the radiation energy 
generation rate is proportional to the local radiation energy 
density. We increased a to 0.1 in this test (only) to reduce the 
thermal timescale, thereby avoiding a long-term drift in the ra- 
diation pressure which occurred only in the 1-D simulations. 1 

The growth timescale predicted by Shakura & Sunyaev 
(1976) in the long- wavelength, radiation-dominated limit is 
5(aQ,)~ 1 = (5/4)t t herm — 8td yn , numerically equal to 2.5 x 10 7 s 
for a — 0.1 and r/r g — 100. Our measured growth time was 
2.5xl0 7 s, a remarkable confirmation of their analytic predic- 
tion. 

On the other hand, in 2-D, convection becomes possible, so 
we now describe Phase B of the simulations. To test for ther- 
mal instability, the disk was first evolved over several thermal 
timescales with the heating rate proportional to the total mass 
(Phase A described in the previous section). We began with 
this prescription for the heating rather than the a-model for 
two reasons. First, we wished to separate evolution driven by 
the convective instability from evolution driven by thermal in- 
stability. Second, as a result of the heat lost by convection, the 
total pressure in the disk after convection sets in is smaller than 
the initial total pressure. Consequently, if the heating rate is 
set by proportionality to the total pressure, the proportionality 
constant that gives thermal balance for the initial disk config- 
uration will not result in enough heating to maintain thermal 
balance after the convective readjustment. Equation [n] ensures 
thermal balance at the beginning of Phase B. 



At the start of Phase B, we first let the thermal instabil- 
ity grow out of numerical noise. Comparing the total radi- 
ation energy in the box at the time the heating prescription 
was changed to a dissipation to the time-averaged radiation 
energy in the simulations with dissipation proportional to gas 
density, the integrated numerical fluctuation in the radiation 
energy was less than 0.1%. In both simulations, as soon as the 
heating prescription was changed, departures from equilibrium 
began to grow. The initial growth was well described by an 
exponential; in simulation 1, its e- folding time was ~ 18td yn , 
~ 0.9 thermal times as measured in the convective equilibrium. 
Simulation 2 showed the same growth time, ~ 18td yn - Because 
the single-zone linear theory of Shakura & Sunyaev (1976) was 
specific to their equilibrium, we cannot directly compare this 
growth time with analytic predictions; however, it is compara- 
ble to the thermal time of the collapsed disk, ~ 16td yn - Thus, 
as predicted by Robertson & Tayler (1981), convection does 
not quench the thermal instability. However, as we shall see 
momentarily, it does significantly change its character. 

In both simulations, numerical noise lead to disk collapse. 
However, the pure a-model without convection predicts that 
the instability has the same growth rate whether the sign of 
the temperature perturbation is positive or negative. To test 
whether radiation energy growth can also occur, we reran Phase 
B with the change that we artificially increased the radiation 
energy density by 5% at the start. In the r = 100r s case, the 
disk still collapsed on a thermal timescale. However, in the 
r = 200r 9 case, the radiation energy density grew on a thermal 
timescale. The simulation was stopped after most of the mass 
was lost through outflow from the simulation grid. We specu- 
late that the positive energy perturbation at r = 100r 9 was re- 
versed because the initial growth in radiation energy triggered 
such strong convection that soon the disk was losing more heat 
than it gained. Because the primary distinction in physical 
conditions between the two cases was that the ratio of radia- 
tion to gas pressure was larger by a factor of 1.7 at r = 200r 9 
than at r — 100r 9 , it is possible that their contrary fates were 
related to this fact. However, we have not been able to con- 
vincingly identify any particular diagnostic that enables us to 
predict which cases will collapse rather than expand despite 
initially positive temperature perturbations. 

We also tried the further experiment of imposing a factor of 
two increase in the local radiation energy density at the end of 
Phase A of the r = 100r g simulation. This had the result of 
the disk heating and expanding until a significant amount of 
mass was lost due to strong outflow. 

On the basis of these numerical experiments, we conclude 
that convection biases the development of the thermal insta- 
bility in the linear regime so that negative radiation energy 
density perturbations are favored. In some circumstances, this 
bias is strong enough that even when the initial perturbation is 
an increase in the radiation energy density, the ultimate result 
of the instability is collapse. However, this is not always possi- 
ble, and sufficiently large positive perturbations can overcome 
convection and lead to runaway growth in the radiation energy 
density. Our simulations, unfortunately, are not capable of fol- 
lowing runaway growth very far because so much of the disk 
mass is quickly pushed out of the problem area. Consequently, 
we cannot say how these cases develop in the long-term. 

In all these numerical experiments we have restricted our- 
selves to very simple phenomenological descriptions of the heat- 



The drift is an artifact of the opacity floor. When the radiation flux is purely diffusive, only the radiation energy gradient matters; conse- 
quently, the general level of radiation energy density is determined only up to an additive constant. In real disks, the diffusion approximation 
breaks down in the optically thin region, so it is possible to think of this additive constant as determined by the condition that the flux must 
be carried by free-streaming radiation at the top of the box. However, in our f-D simulation the flux is diffusive everywhere due to the opacity 
floor, so there is nothing in the equations to prevent numerical drifts in the general level of radiation energy density. By contrast, in the 2-D 
simulations, the flux leaving the box is carried by advection of radiation, which then fixes the radiation normalization as in the free-streaming 
case. 
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ing and, in the case of the simulations with imposed pertur- 
bations, very simple and uniform perturbations. Real disks 
undoubtedly behave differently. We regard these simulations 
therefore as demonstrating that it is possible for thermal insta- 
bility to occur after convective readjustment when the heating 
rate is roughly described by the Q-model, and that in some cir- 
cumstances, convection can cause a bias toward collapse. How- 
ever, these simulations certainly do not provide the final answer 
to the question of whether and how thermal instability affects 
physical disks. 



4.2. Final state 

After the initial exponential collapse, the radiation energy 
density approaches a constant value as the disk becomes gas- 
pressure supported and thus thermally stable. This is a con- 
sequence of the fact that for fixed surface mass density E and 
fixed stress-parameter a, there exist exactly two thin disk so- 
lutions, one radiation-pressure supported with high accretion 
rate and flux, and one gas-pressure supported with low accre- 
tion rate and flux (an advection-dominated solution can also 
exist, but it is no longer thin). Figure ^ shows the radially 
averaged gas density in the final state of simulation 1, phase B, 
while figure (7] shows the radiation and gas pressures. Despite 
the fact that the gas pressure is is larger within the bulk of the 
disk, Prad becomes larger than p ga3 near the photosphere sim- 
ply because the gas density drops precipitously. The density 
weighted mean of \z\ is smaller by a factor of 12.5 at the end 
of Phase B than at the beginning of Phase A. Twice the grid is 
rebinned by a factor of two to follow the collapse. We find that 
the gas becomes strongly clumped in the radial direction filling 
only ~ 1/4 — 1/2 of the simulation region. Since centrifugal 
forces balance gravity in the radial direction, the clumps re- 
main stable. In a real disk, we expect that viscous stresses will 
erase this radial stratification. 

To illustrate the relationship between the initial radiation- 
pressure supported state and the final gas-pressure supported 
state, we explicitly find them in terms of the parameter (j> = 
Prad/Pgas = (|aT 4 ) / ' (pk bT / ' m) (m is the average mass per 
particle). Ignoring the radial clumping, we approximate the 
disk as a single vertical zone from z = to H , writing the 
hydrostatic equilibrium equation in the form: 



(1 + 4>)Pgas = pH 2 fl 2 . 



(22) 



In the q prescription, the radiation equilibrium equation be- 
comes 



3 , J.N tr 4:C4>Pgas 

~afi(l + 4>)pgasH = 

2 acL 



Defining E = 2pH, we solve for 



(l + ^) io = 6 7o 



where 7 is defined as 

la -( 8c v 

\3QK es E/ 



2am 



(23) 



(24) 



(25) 



When radiation pressure dominates, cf> ^> 1, so 4> ra d ~ 7o 



1/4 

O 

-1/6 



When gas pressure dominates, if </> 0.1, then 4> ga s ~ 7o 
Thus, we expect that 4> gas ~ 4>~^ Z ■ If 4' aas ~ 1) then it must 
be solved for numerically using equation E4. 
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Fig. 6. — Radially-averaged density profile at the end of simulation 
1, phase B. 
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Fig. 7. — Radially-averaged pressure profiles at the end of simu- 
lation 1, phase B. Solid line is radiation pressure, while dashed line 
is gas pressure. 

Due to the clumping, this relation is only qualitatively borne 
out by our simulations. Simulation 1 begins with <j> ra d = 95 
(density- weighted average). As a result of the a renormaliza- 
tion that occurs in order to preserve thermal balance when the 
heating prescription is changed, a changes from 0.01 to ~ 0.05 
(since the radiation pressure is reduced by ~ 1/5). With this 
new value of a, we would expect on the basis of the simple 
one-zone model just presented that the final 4>gas after collapse 
should be ~ 0.9; we find 0.54 in simulation 1. Similarly, sim- 
ulation 2 begins with <j> = 160. The effective a again changes 
from 0.01 to 0.05, leading to a predicted final 4> ga s = 0.38 
(from equation |2J); we find 4> ga s — 0.95. The column density 
is artificially increased due to the clumping, so the effective 
7o oc E -8 decreases (equation ^BJ), leading to a larger 4>gas- 
The radial clumping is much stronger in simulation 2 than in 
simulation 1, which explains the larger discrepancy. In simula- 
tion 2, p r ad ~ Pgaa in the final state, so weak convection is still 
present. 

5. CONCLUSIONS 

5.1. Evolution on the dynamical time- scale 

We have shown that in a very short time (of order the dy- 
namical time), convection modifies the structure of geometri- 
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cally thin radiation-dominated disks. Although there are or- 
der unity contrasts in specific entropy at fixed altitude z, the 
radially-averaged entropy becomes very nearly constant with 
height as a result of convective mixing. The mean density pro- 
file that emerges is close to the one predicted in the analytic so- 
lution of BKB. As likewise predicted by that analytic solution, 
the convective flux near the disk midplane is comparable to the 
diffusive flux. More rapid upward heat flux (for the same dissi- 
pation rate) results in smaller mean radiation pressure, weaker 
support against gravity, and therefore greater mean density. 
The decreased pressure may cause a further decrease in the 
viscous stress, resulting in a higher surface mass density and 
thus lower radiation to gas pressure ratio. One consequence 
may be that the region of thermal instability is reduced in size 
by a factor of a few, changing the accretion rate at which near- 
Eddington disks are subject to thermal instability. 

5.2. Evolution on the thermal timescale 

We have also shown that evolution on the thermal timescale 
is quite sensitive to the specific character of the dissipation pre- 
scription. Dissipation proportional to density leaves the disk 
in a long-term statistical equilibrium in which convection and 
radiative diffusion carry almost equal amounts of energy. 

On the other hand, dissipation proportional to vertically- 
integrated pressure is subject to the thermal instability first 
pointed out by Shakura & Sunyaev (f976), but subject to a 
significant modification due to convection: it is biased toward 
collapse, rather than runaway expansion. 

The result of thermal collapse is a new gas pressure- 
dominated equilibrium at the old column density. The pres- 
sure in this new equilibrium is much smaller than in the orig- 
inal one. Because the stress such a disk is capable of exerting 
(assuming the a model still holds at least approximately, as ex- 
tensive MHD simulations indicate: Balbus & Hawley f998) is 
much smaller than in the initial state, the mass accretion rate 
through such a disk is much smaller than the initial value. Con- 
sequently, if outer regions of the disk continue to pass matter 
inward at the same rate, mass must build up in the collapsed 
portions of the disk. 

The long-term evolution (i.e., on the inflow timescale) of this 
situation is difficult to predict. We speculate that the column 
density will continue to build up until thermal instability sets in 
at Prad ~ Pgas- The specific point at which this occurs may be 
modified by convective heat transport. At this point, the disk 
no longer has a gas-pressure supported equilibrium available 
for fixed column density, so the disk may heat up, joining the 
advection-dominated "slim disk" solution, leading to limit cy- 
cle behavior on the inflow timescale (Abramowicz et al. 1988), 
or convection may lead to episodic release of energy on the dy- 
namical timescale, matching on average the energy dissipated 
within. 



5.3. Future improvements 

We have made several assumptions to simplify our calcula- 
tions that can be tested with future work. First, we have as- 
sumed that flux-limited diffusion is an appropriate description 
of the radiation field. Since our simulations were optically thick 
everywhere, this may be appropriate, but needs to be checked 
with a more accurate computation, for example using the vari- 
able Eddington factor method described in Stone, Mihalas, & 
Norman (1992). 

Second, we have assumed simple dependences of the heating 
rate on local disk conditions. Actual disks may have a heating 
rate which depends on global disk parameters, e.g. height or 
radius. If a large fraction of the magnetic energy is carried to 
the corona, then the disk may never become radiation-pressure 
supported (e.g., as suggested by Svensson & Zdziarski 1994). 
Simulations of disk annuli indicate that a significant fraction 
of the magnetic energy generated within an accretion disk can 
be carried to the corona (Miller & Stone 2000), but not neces- 
sarily enough to eliminate a radiation-pressure dominated disk. 
These simulations need to be extended to include radiation and 
a physical equation of state to determine the magnitude of the 
transported energy. Because disks with dissipation primarily 
at the surface are gas pressure-supported, they are subject to 
neither convective nor thermal instability. 

We believe that the problems of the previous paragraph may 
best be addressed with 3-D radiation magnetohydrodynamic 
simulations in which the magnetic field strength and dissipa- 
tion rate are computed self-consistently. The nature of con- 
vective transport as well as the contribution of turbulence to 
heat transport (Balbus 2000) may have a large effect on our 
results. Magnetohydrodynamic turbulence might destroy the 
convective plume structure since the magnetorotational insta- 
bility operates on the dynamical timescale. The damping of 
turbulence by radiation diffusion or shocks at the scale aH 
might determine the heating profile of the disk (Agol & Krolik 
1998). If the magnetic stress scales with gas pressure rather 
than radiation pressure, the thermal instability is suppressed 
(Piran 1978); whether this is the case can be addressed with ra- 
diation MHD simulations. We have purposely limited the scope 
and detail of our current simulations as our neglect of 3-D mag- 
netohydrodynamics will change the nature of the equilibrium. 
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